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Abstract. Single Helical Axis (SHAx) configurations are emerging as the 
natural state for high current reversed field pinch (RFP) plasmas. These states 
feature the presence of transport barriers in the core plasma. Here we present 
a method for computing the equilibrium magnetic surfaces for these states in 
the force-free approximation, which has been implemented in the SHEq code. 
The method is based on the superposition of a zeroth order axisymmetric 
equilibrium and of a first order helical perturbation computed according to 
Newcomb's equation supplemented with edge magnetic field measurements. The 
mapping of the measured electron temperature profiles, soft X-ray emission 
and interferometric density measurements on the computed magnetic surfaces 
demonstrates the quality of the equilibrium reconstruction. The procedure for 
computing flux surface averages is illustrated, and applied to the evaluation of 
the thermal conductivity profile. The consistency of the evaluated equilibria with 
Ohm's law is also discussed. 
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1. Introduction 

The reversed field pinch (RFP) configuration has been known for a long time as a 
possible candidate for the magnetic confinement of fusion plasmas [3 [2]. Among 
toroidal configurations, it is the only one which could in principle achieve reactor- 
relevant conditions without the need of additional heating systems, thanks to its 
high level of ohmic power dissipation. Nevertheless, up to recent years the RFP was 
regarded as having poor confinement properties, due to the simultaneous presence 
of several MHD tearing modes with m=l poloidal periodicity 3 . These modes are 
resonant at different radii within the plasma core, and their magnetic islands overlap, 
so that magnetic surfaces are destroyed and a high level of magnetic chaos ensues. As a 
consequence, plasmas with flat density and temperature profiles over the whole plasma 
core are produced in this condition, which is dubbed Multiple Helicity (MH) state. 
The presence of these modes is intrinsic to the configuration, since they are required 
to drive poloidal currents in the outer plasma region while the applied electric field 
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is only toroidal (dynamo effect) [H IS]- In the past, chaos reduction was obtained 
transiently by inductive poloidal current drive [5J [7] , but a way of obtaining this effect 
in a stationary fashion is still lacking. 

This view of the RFP as a stochastic plasma with bad confinement properties is 
radically changing, thanks to the recent discovery of a new class of equilibria, which 
have been dubbed Single Helical Axis (SHAx) states [HI [9]. SHAx states are the result 
of a natural tendency of RFP plasmas to move, as dissipation changes, towards a 
single helicity (SH) condition, that is a condition where only one of the core-resonant 
tearing modes provides the dynamo effect (TUl El dH US] • This tendency, which has 
been predicted since several years by 3D MHD simulations [M], has been confirmed 
by experimental results in several RFP devices [151 HB HZl [H] , where states in which 
the innermost resonant m=l mode dominates over the others (secondary modes) have 
been observed. Such conditions are generally called Quasi Single Helicity (QSH) states, 
since they differ from the theoretical SH ones due to the fact that the secondary modes 
still have a non-negligible amplitude. In particular, in the RFX-mod device (R = 2 
m, a = 0.459 m) it has been observed, when feedback control of the radial magnetic 
field at the edge is applied [lH [20] , that the duration and frequency of occurrence of 
QSH states increases with plasma current [9l[2TJ[22]. The secondary mode amplitude, 
normalized to the total magnetic field, is also found to decrease with increasing plasma 
current, so that experimental QSH states are progressively approaching the ideal SH 
condition [9l [21] [22] . Plasma current increase is strongly correlated to an increase of 
electron temperature |23j . 

In the QSH states the dominant mode amplitude is found to increase with plasma 
current, inducing a bifurcation to a new topological structure of the magnetic field 
[S]. This happens because, beyond a threshold in the dominant mode amplitude, the 
X-point of the magnetic island separatrix collapses on the main magnetic axis, and the 
two disappear, leaving the original island 0-point as the only magnetic axis. Thus, 
a helical plasma column is obtained in the axisymmetric device. These SHAx states, 
which represent a special flavour of the more general QSH condition, are more resilient 
to magnetic chaos, due to the disappearance of the separatrix They constitute a 
new paradigm for high performance RFP plasmas, which could lead to a re-evaluation 
of the potential of this configuration for fusion reactor development. 

In this paper we propose a method to determine the flux surface shape of RFP 
discharges in a SHAx state, which involves only the solution of ordinary differential 
equations. The method, which has been implemented in a code named SHEq (Single 
Helicity Equilibrium) starts from the symplectic representation of the magnetic field 



where (p, 9, (p) is a generic system of curvilinear coordinates used to describe the 
toroidal system. F = F{p, 6, cf)) and = '^{p, 0, (j)) represent the toroidal flux through 
a constant F surface and the poloidal flux outside a constant 5* surface, respectively, 
both divided by 2n. The existence of magnetic flux surfaces is assured if a function 
jo(x) such that 



exists. Flux surfaces are then described by the equation p = const, (or f{p) = const. 
for any / function of p only). The magnetic axis is defined by V p = 0. In the case 
of an axisymmetric plasma with circular cross-section the function p can be taken to 




B = Vi^ X V6I - V^- X V(/> 



(1) 



B • Vp = 



(2) 
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coincide with the radius r of the magnetic flux surfaces, and F(r) and ^'(r) may be 
interpreted as the usual toroidal and poloidal flux, respectively. In the case of a general 
magnetic field, F{r, 9, (f) and 4'(r, 9, (p) in (1) are not constant over flux surfaces, flux 
surface existence being not guaranteed in the general 3D case. In other terms, we can 
say that adding a generic perturbation field to the axisymmetric one, the circular flux 
surfaces are destroyed, and it is not clear a priori if some other flux surfaces exist, and 
which can be the function p that labels them. 

If one assumes a helical symmetry of the perturbation, so that and F depend 
only on r and on m = m9 — ncf), the symmetry ensures flux surface existence, and it 
is straightforward to show that the helical flux x(r, u) — nVl!{r, u) — nF(r, u) is a flux 
function, that is B • Vx — 0. Thus, the knowledge of x('';'") enables to display the 
shape of the magnetic surfaces. This is the basic hypothesis adopted to treat SHAx 
states. 

The method that we propose for computing the helical flux in SHAx states is based 
on the approach described in ref. |26] for calculating the tearing mode eigenfunctions 
in a force-free RFP plasma in toroidal geometry. The SHAx states are therefore 
considered as composed of a dominant, saturated tearing mode superposed on an 
axisymmetric equilibrium. The tearing mode amplitude is assumed to be small with 
respect to the axisymmetric fields, so that a perturbative approach can be adopted. 
The same coordinate systems and the same notation of ref. [55] will be used. The 
most relevant results of ref. [26] are reviewed in the two sections next, so as to give 
a self-consistent treatment of the equilibrium problem, the reader is referred to that 
paper for further details. 

This paper is organized as follows: in Section [2] the zeroth-order axisymmetric 
equilibrium calculation is introduced; in Section|3]the first order correction due to the 
dominant mode is computed; in Section|4]the method used to compute the flux surfaces 
is described and validated by studying the mapping of electron temperature, density 
soft X-ray emission measurements on them; in Section [5] a method for computing 
flux surface averages is described, and an example application to power balance 
computation in a SHAx condition is given; Section [S] is devoted to a discussion of the 
constraint given by the flux surface averaged Ohm's law parallel component. Finally, 
in Section [7] conclusions are drawn. 

2. Zeroth-order equilibrium 

Let us consider a zeroth-order axisymmetric toroidal plasma with circular cross- 
section, formed in a vacuum chamber with major radius Rq and minor radius a. The 
flux surfaces are non-concentric circles, each having radius r, being horizontally shifted 
by a quantity A(r). The shift of the outermost flux surface is imposed as boundary 
condition (in the experiments this is obtained from external magnetic measurements). 
A point lying on one of these flux surfaces is identified by the radius r of the surface, 
by the poloidal angle 9 measured with respect to the inboard mid plane, and by the 
toroidal angle <j). These coordinates, which we call geometric coordinates, u*, are 
related to the standard cylindrical system (i?, Z) used to describe toroidal fusion 
devices by 



i? = i?o-?'cose' + A(r) 
Z = rsin^. 



(3) 
(4) 
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where Rq is the torus major radius. The (r, 9, 0) coordinate system is curvihnear 
and non-orthogonal, in order to properly take into account the toroidal geometry. 
A complete knowledge of the metric tensor is essential (see Appendix A for a brief 
reminder on curvilinear coordinates). The metric tensor and the Jacobian of the 
geometric coordinates are given in Appendix B. 

The contravariant representation of the zeroth-order magnetic field associated 
with the geometric coordinate system is 

Bo = VFo(r) xV0- V^olf) x + Vr x Viy{r, 6) (5) 

where Fq and are, respectively, the toroidal and poloidal flux divided by 2tt. The 
equilibrium is fully defined once -Fo(^)j ^o('')i ^('') ^{rjO) are known. 

Following the standard procedure for introducing fiux coordinates [27j . one can 
define a new poloidal angle as 

d = e + X{r,9) (6) 

with X{r,9) = v{r^9) / FQ[r). Here and in the following a prime designates derivative 
with respect to r of quantities which are functions of r only. In the = (r, i?, (f) 
system, which we call flux coordinates, the magnetic field lines are straight and the 
magnetic field contravariant representation is simply 

Bo = VFo(r) X Vi?- V*o(0 X V0. (7) 
This provides simple formulas for the contravariant components -B', given by 

55 = Bt = ^^', Bt^-^FI, (8) 

where ^/g^ is the Jacobian of the flux coordinate system: 

^ = (Vr • V?9 X V0)"\ (9) 

The metric tensor and the Jacobian of the flux coordinates are also given in Appendix 
B. 

The determination of the parameter X{r,9) in ^ is possible for a large aspect 
ratio torus following a perturbative approach. Ampere's law allows to deduce the 
current density contravariant components J'. From the force balance condition and 
Bq = one gets Jq = 0. Using this information, and performing an expansion in the 
small aspect ratio parameter e = a/i?o, it is possible to compute the quantity relating 
9 and as 



where 



X{r,9) = Ai(r)sine' + A2(r)sin26' + o(e^). (10) 

Ai(r) = ^-A'(r) A2(r) = -^Ai(r). (11) 
The inverse of transformation ^ is then easily derived as 

6l = z?-Aisini?- (^A2 - sin2i9 + o(e3). (12) 

Using equation (ITU|) the relation between cylindrical {R, (f), Z) and flux coordinates 
can also be found; 

i? = i?o - r cos i9 + A(r) - rXi (r) sin^ i9 + f A^ - 2r A2 j sin^ cos + o{e^h) (13) 
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Z = rsmiS- ^Ai(r)sin2i9+ i^^rXf - 2r\2^ sin^Jcos^ - -A^(r) sint? + o(e^&) (14) 

so i? EE R{r. {)) and Z = Z{r, d). 

The flux coordinate Jacobian for a large aspect ratio torus is, as stated in 
Appendix B, 



with 



(15) 



K{r) = -^{l + ^ + -^A' -^ + o{e^)) . (16) 



^ [ 1 + A + - — ^ r^frS- 

r \ Ro 2Ro 2Rl 

It is possible to show that, for a force-free equilibrium, in flux coordinates the current 
density is proportional to the magnetic field through a coefficient which is a function 
of r only, that is 

MoJo = CT(r)Bo. (17) 

It is convenient to define, for a generic field A, its hatted version as A = ^/g^A, 
which hides the Jacobian contribution. The zeroth-order hatted magnetic field and 
current density components, function of r only, are then: 

B^ = % Bt=F^ f^oJ^ = a^'o f^oJ^ = cTF^. (18) 
Given the cr(r) profile, which is an input to the algorithm, the zeroth-order force 
balance yields the following equations: 

^[K{r)Bt] = -air)Bt (19) 



dr 
d_ 
dr 



-Bq 



B^=a{r)Bl (20) 



The second equation contains metric coefficients which are function of r and By 
using the expansion in harmonics described in Appendix B, it can be split into two 
equations, one for Bq and one for A. Furthermore, because of the nonlinearity given 
by the fact that the metric coefficients depend on A, it is convenient to introduce a 
perturbative expansion: 

B^^Bf + B^..., Bt = o{e^)Bt (21) 

Bt = Bt + Bt + .... Bt = o(^)Bt. (22) 

The resulting equations for the lowest order contribution are 

d f Ro 



d_ 

Jr \R 



' Br]=air)Bf (24) 



which can be solved for Bf(r) and Bf{r). The solution starts from the magnetic axis, 
where regularity imposes Bf ~ aoRor/2 and Bf ~ r, and proceeds to the edge. The 
solutions are then rescaled so as to match a boundary condition, for example the edge 
poloidal field. 

These solutions are then be plugged into the equation for the shift 
, „ dBf/dr\ 

l + 2r— 1^ 1 +^=0. (25) 
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which is solved using the boundary conditions A'(0) — and an assigned value for 
A(6), b being the radius at which the magnetic measurements yielding the condition 
are located. 

Having determined the shift A(r), the next order correction to the fields can be 
computed by 



d_ 
dr 
d_ 
dr 



Rq 



d_ 

dr 
d 
dr 



A 



r 

i?0 



2Rl 



r 

2Ro' 
~2~ 



2Ro -Ro 



(26) 



(27) 



Again, these equations are solved starting from the axis, where regularity requires 
~ 3r/(2CToi?o) and B^ ~ [3/(a^R^) - Ao/Ro]r, Aq being the shift of the magnetic 
axis. In practice this correction turns out to be very small. 

It is worth noting that this method of computing the zeroth-order axisymmetric 
equilibrium, if compared to the standard Grad-Shafranov equation, has the advantage 
of requiring simply the solution of five ordinary differential equations. This is obtained 
at the price of being restricted to deal with circular flux surfaces, which is however 
reasonable for present day RFP devices. While cr(r) can in principle be any function, 
for the application described in the following the customary parametrization called 
a-Qo model has been used, that is 
260 



1 



(D 



(28) 



The two free parameters Qq and a are adjusted so as to obtain given values of the 
two dimensionless parameters Q and F, which are the well known pinch and reversal 
parameters used to describe RFP plasmas. 



3. First-order contribution to equilibrium 

The next step is to add a non- axisymmetric perturbation to In the gauge Ar = 0, 
where Ar is the covariant radial component of the vector potential A, the total 
magnetic field can be written as 

B = VF X V?9 - V^- X V(/> (29) 

where now F and ^> depend on all three coordinates. The = (r, i?, 0) coordinates 
are not any more flux coordinates for the perturbed magnetic field, which means that 
magnetic field lines of B are not straight in this coordinate system, and the potentials 
F and 4" are not any more flux functions. These potentials, related to the vector 
potential covariant components, can be Fourier expanded as 

vE'(r,^,^) = -^^(r,^,0) = *o(r)+ ^ ^™.«(^y(m^-n0) (gg) 

F(r,^,0) = ^4r,,9,0) = Fo(r)+ ^ (31) 

The perturbed quantities contain n 7^ terms only, and the harmonics amplitudes are 
complex. 

Given the representation l|29|) of the magnetic field, the total hatted contravariant 
magnetic field components are 

^,0 9^- dF dF d"^ 
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Computing the total current density components from Ampere's law and plugging 
them into the first-order force balance equation 

j X Bo + Jo X b = -^e'^'^CfH + J^V^w^ = (33) 



9n 

one obtains the proportionality between perturbed radial current and perturbed radial 
magnetic field 



IJ-of = cr(r)6'^ 



and 



d 



d 



d-d ^d<j) 



b^^^O 
dr 

b q— = 0. 

dr 



(34) 

(35) 
(36) 



By Fourier-transforming equations (1341) and p6l) (only two equations are needed, 
since for each mode there are two unknown functions f™--" and V'™'") and using 
Ampere's law the following equations are found: 



dr 



0,0 



~in 



9 m 



[nf 
'd%l;' 



dr 



,m+l,n 



dr 



nf 

dtp'' 



(m + 1)^ 



(m — l)ip 



771 — 1,711 



dr 



(37) 



/9n ' 



-m 



dtp' 



771 + 1,71 



dr 



dr 



rij ' — mip — ■ 



/gw J ' ' m — nq dr 

(38) 



The method used for the solution of these equations is described in detail in ref. 
|26| . The solution requires the knowledge of the corresponding harmonic amplitude 
for the radial component of the magnetic field at some surface outside the plasma, 
which represents the boundary condition. Furthermore, if the mode has a resonant 
surface inside the plasma, a discontinuity in the eigenfunction derivative should be 
allowed. The magnitude of this discontinuity is obtained by imposing a further 
boundary condition, that is the amplitude of the toroidal magnetic field component 
at the same surface where the radial one is determined (that is the surface where the 
sensors are located). 
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4. Helical flux surfaces 

As already mentioned in the introduction, the SHAx states will be modeled as pure 
Single Helicity states composed of the superposition of the zeroth-order axisymmetric 
equilibrium and of the dominant mode eigenfunction. The dominant mode, in the 
case of RFX-mod, is the m=l/n=7. An example of the m = 1 mode spectrum for a 
1.5 MA SHAx state obtained in RFX-mod is shown in Fig. [TJ The spectrum gives the 
amplitude of the toroidal field component, normalized to the average poloidal field, 
measured outside the plasma. It can be seen how the n = 7 mode clearly dominates 
over all the others, justifying the choice of taking it to be part of the equilibrium. 




Figure 1. m = 1 mode spectrum in a 1.5 MA SHAx state in RFX-mod (sliot 
24598, t = 174 ms). Tiie amplitudes of tiie Fourier modes on ttie toroidal field 
component, normalized to the average poloidal field, are shown. 

The poloidal and toroidal flux will be taken as: 

*(r, u) = ^Q{r) + V'"'"(r)e™ + c.c. (39) 
F(r, u) = Fo{r) + /'"■"(r)e™ + c.c. (40) 

where 

u = md — ncj) (41) 

and ^'0(7'), Fo{f), -(/^"'"(r) and f'^^'^ir) are computed as described in the preceding 
sections. The helical flux 

x(r, u) = m^E' — nF ~ 

= mfo - nFo + (mV'"'" - n/"^") exp(m) + c.c. (42) 

is constant on the resulting flux surfaces. An example of such surfaces, obtained as 
contours of x, is shown in Fig. [2] It can be seen that only the inner surfaces are 
significantly distorted, and assume a bean-like shape. On the contrary, the outer 
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Figure 2. Flux surfaces in a 1.5 MA SHAx state in RFX-mod (shot 24598, t = 
174 ms), computed as contours of the helical flux x- 



ones retain a quasi-circular shape, with a shift induced by the presence of an m = 1 
component in the equihbrium. 

SHAx states often display profiles of the kinetic plasma quantities which are non- 
symmetric with respect to the vacuum chamber axis. This is a direct consequence of 
the asymmetry of the helical flux surfaces displayed in Fig. [5] A direct test of the 
goodness of the equilibrium reconstruction is therefore a study of the profiles of such 
quantities when plotted as a function of a flux function. In the following we shall use 
as abscissa for such plots the effective radius defined as 



^ — ^ (43) 

Xa - XO 

where xo is the helical flux on the discharge helical axis and Xa is the helical flux of 
the outermost surface. The p variable ranges between (on the helical axis) and 1. 

The first test of our equilibrium reconstruction is performed by considering an 
electron temperature (Tg) profile measured along a horizontal diameter of the chamber 
by an 84-point Thomson scattering system [28j . An example of the measurements is 
shown in Fig. [3^, where the previously mentioned asymmetry can be clearly seen. It 
is also observed, as already reported elsewhere jJI], that the profile features strong 
electron transport barriers, revealed by the steep gradients in the shaded regions of 
the graph. The same measurements, plotted as a function of p, are shown in Fig. 
[Sjj. It is immediately apparent how the two half profiles, which in the previous plot 
were evidenced by the use of blue and red points, collapse one onto the other. This 
is a proof that the helical flux evaluated with our method is indeed a flux function, 
assuming that Tg is (as reasonable given the very fast parallel thermal transport and 
the fact that the density proflle is essentially flat). Fig. [5J: shows the resulting electron 
temperature map on the poloidal plane. The hot bean-shaped central region is clearly 
seen. 
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Figure 3. (a) Electron temperature measured along a horizontal diameter of 
RFX-mod by the Thomson scattering technique. The two colours mark points 
that are on the two sides of the geometric chamber axis, (b) The same data of 
frame (a) plotted as a function of normalized the helical flux p. It can be seen 
how the two halves of the temperature profile collapse one onto the other, (c) 
Temperature map on the poloidal plane. The data refer to shot 24599 at t = 99 

ms- (10"Wm") 

A second test has been performed using a 78-chord tomographic system, which 
measures the soft X-ray emission of the plasma In this case, being the 

measurements integrated along hnes of view, the emissivity e has been assumed to 
be a function of p only through a simple four-parameter model: 

e(p) = eo(l - p")^ + £1 (44) 

The four parameters have been optimized so as to minimize the sum of the squared 
differences between the actual measurements and those reconstructed from the model. 
The final emissivity profile is shown in Fig. H^. The two classes of data, measured and 
reconstructed, are reproduced in Fig. |4)d. The figure displays two sets of points for 
each class, corresponding to the three vertical fans and to the single horizontal fan of 
measurement lines, respectively. It can be seen that the match is very good, despite 
the simplicity of the used model. The resulting emissivity pattern in the poloidal 
plane is shown in Fig. |4t, showing once again the hot bean-shaped core region. 

(10"Wm") 




0.2 0.4 0.6 0.8 1.0 -0.4 -0.2 0.0 0.2 0.4 -0.4 -0.2 0.0 0.2 0.4 

P r (m) r (m) 



Figure 4. (a) Soft X-ray emissivity profile given by the model of Equation 1144 l l. 
after parameter adjustment, (b) Measured values of line-integrated soft X-ray 
emission (red diamonds) and values reconstructed according to the emissivity 
profile shown in panel a (black triangles). Emissivity pattern on the poloidal 
piano (c). The data refer to shot 26608 at t = 160 ms. 

Finally, line-averaged density measurements given by a 7-chord two-color 
interferometer [30] have been analyzed. Once again, a simple model, in this case 
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with four parameters, has been assumed for the density dependence on p: 

n{p) = no - (no - ni - nedge)p°' - nip'^ (45) 

with Uedge = 10^* m~'^. Usually, density profiles are flat, due to the fact that the 
discharge fueling is done by hydrogen released by the graphite first wall or puffed with 
valves. Thus, the density profile analysis has been performed on a case where a pellet 
had been injected into the plasma. The density profile resulting from the optimization 
of the model parameters is shown in Fig. [5^, showing a hollow profile due to the efl^ect 
of pellet ablation in the outer portion of the plasma. This is is different than the 
example shown in Ref. !9!, where a peaked profile resulted from the pellet ablation 
taking place in the inner bean-shaped region. The measured values of line-integrated 
density are shown in [SJd, together with the reconstructed values. A good agreement 
can be seen also in this case. Finally, the density pattern on the poloidal plane is 
reproduced in [5];. 



(10" m'^) 




0.0 0.2 0.4 0.6 0.8 1.0 -0.4 -0.2 0.0 0.2 0.4 -0.4 -0.2 0.0 0.2 0.4 

P r (m) r (m) 



Figure 5. (a) Plasma density profile given by the model of Equation 1)45 p . after 
parameter adjustment. (b) Measured values of line-integrated density obtained 
through microwave refiectometry when a hydrogen pellet enters the bean-shaped 
plasma region (red diamonds) and values reconstructed according to the density 
profile shown in panel b (black triangles). Density pattern on the poloidal plans 
(c). The data refer to shot 24936 at t = 88 ms. 

Summarizing, we have thus demonstrated, theoretically and experimentally, that 
the helical flux is a good fiux function for SHAx states. Furthermore, the good 
mapping of the kinetic quantities on the computed fiux surfaces shows that such 
surfaces are indeed, at least partially, found into the plasma. This is a proof of the 
emergence of some degree of order from the magnetic chaos present in MH condition. 

5. Flux surface averaging 

The computation of flux surface averages of different quantities is an essential step 
for building transport equations depending on one coordinate only (the flux surface 
label). The coordinate systems used up to this point are not appropriate for describing 
the helical fiux surfaces, due to the fact that there are points which are not uni vocally 
identified by them. This can be understood by considering that for any inner fiux 
surface, which does not contain the magnetic axis of the axisymmetric equilibrium, it 
is possible to identify two points which have the same value of the poloidal angles 6 
and 'd. We thus need to build a new coordinate system. This is done by adopting x as 
radial coordinate, keeping (j) as toroidal angle, and introducing a new poloidal angle. 
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f3, which rotates around the hehcal axis (Vx = 0). The definition of /3, with respect 
to the cyhndrical coordinates, is 

where Ra{4>) and Za{(t>) represent the coordinates of the hehcal magnetic axis. We refer 
to Appendix C for the derivation of the metric tensor of the new /3, (p) coordinates. 

The flux surface average of a generic quantity ^(x,/3, in a toroidal system is 
given by |27] 

ffdpd^^g- ^^'^ 

The Jacobian ^ of the (x, /3, 0) system can be computed relatively easily in terms of 
the Jacobian of the zeroth-order flux coordinates w"^ used above. One obtains 



The derivatives appearing in expression (j48p can be computed from the x values given 
by Newcomb's equation and from recalling (IT51) and (ITi|) . 

As an example of the application of this method, we show in Fig. [6] the flux 
surface averages of the toroidal and poloidal covariant magnetic field and current 
density components. Using the relations that link the helical coordinates to the 
system of zeroth-order flux coordinates, the held components in the helical system can 
be obtained from those computed in the u;* system. Equations ([5^ return the hatted 
contravariant component of the magnetic held, so the metric contribution must be 
reintroduced: despite the assumption of single helicity in the fluxes, the contribution 
of the Jacobian harmonics (formula (jB.2ip ) generates non-zero field components also 
for different mode numbers. This happens because of the toroidicity, which produces 
a toroidal coupling between modes with same toroidal mode number n and different 
poloidal mode number m. The components result: 

B'' = i (n/i'"(r) - mV'^'"(r)) G(r, 0) + c.c. (49) 
B^ = f — j— + —1= cosd] %{r) + V''^^"(r)G(r, i?, 4>) + c.c. (50) 

1 2 \ 

: cos i9 F^(r) -K/'i'"(r)G(r, 1^,0) c.c. (51) 



with 

G{r,'d,(j)) — gt(mi?-n0) ^ 1 gt[(m-l)i?-rt0] ^ ^ ^i[{m+l)'d-n4>] ^^2) 

^/9wa \Jgw\ \Jgw\ 

The definitions of g^o and g^x are given in Appendix B. 
Writing 

B = g''er + g'^e.g + g'^e0 , (53) 

Bpoi Bfor 

the magnitudes of the poloidal and toroidal magnetic field components are: 

Btor ^B^ = ^fg^B^ (54) 



Bpol = ^grr{B-y + 2gr^B-B^ + gMB-"?. (55) 
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Btor can also be identified with the measured covariant components on the 
normahzed basis vectors e^. 

The graphs of Fig. [5] show that the flux surface average of the toroidal field 
component is monotonically decreasing, as for the standard cylindrical models of the 
RFP fields, with a maximum of 1.1 T which is now located on the helical axis, a slow 
decrease in the central part of the plasma, a knee around p = 0.25, and a reversal in 
the outer part of the plasma. The poloidal component is also rather flat on the inner 
surfaces, where it has an almost uniform value around 0.5 T. The same features are 
displayed by the current density components, not surprisingly since this is a force-free 
equilibrium. The maximum toroidal current density, on the helical axis, is around 6 
MA/m^, while in the inner part of the plasma the poloidal component takes values a 
little larger than 2.5 MA/m^. It can be also remarked that for p > 0.8 the toroidal 
current density is negligible. 




Figure 6. Left: Flux surface averages of the toroidal and poloidal magnetic field 
components. Right: Flux surface averages of the toroidal and poloidal current 
density components. The data refer to shot 24598, at i = 174 ms. 

Using the same approach to flux surface averaging, one can compute the average 
ohmic power, (77 J^). In doing this the Spitzer-Harm resistivity formula has been 
used, with the electron temperature profile measured by the Thomson scattering 
system. Furthermore, we have assumed a uniform effective charge profile, and the 
correction factor with respect to the Z = I resistivity has been computed comparing 
the volume integral of the dissipated power to the actual input power, P — VI, 
assuming stationary conditions. The result is shown in the left panel of Fig. [T] For 
this case, the correction factor turned out to be equal to 1.75, which appears a rather 
reasonable value. It is interesting to notice how the dissipated power displays a peak 
outside of the bean-shaped inner region. It should however be emphasized that this 
feature crucially depends on the profile chosen for cr(r). 

The averaged dissipated power can be used in a simplified power balance equation, 

to compute the thermal conductivity k. Here g^^ = Vx • Vx is the first metric tensor 
element, which has been computed explicitly in Appendix C, while V' ~ dV/dx is the 
specific volume. Using the temperature profile given in Fig. [3] and the averaged power 
of Fig. [7] (left), the thermal conductivity profile shown in the right panel of Fig. [7] 
is obtained. It can be seen that the thermal conductivity displays a minimum, at a 
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P P 

Figure 7. Left: Flux surface average of dissipated ohmic power. Right: Thermal 
conductivity profile computed from the surface-averaged power balance. The data 
refer to shot 22182, at i = 49 ms. 



value around 8 m^/s, corresponding to the strong gradient in the temperature profile. 
This is an order of magnitude lower than values obtained in MH conditions [3T] . 

6. Ohmic constraint 

As in the RFP there is a large current flowing in the plasma, a natural question 
which arises when an ideal equilibrium has been computed is whether such current is 
consistent with Ohm's law. In MH states this cannot be easily assessed, since Ohm's 
law includes a "dynamo electric field" which results from the effect of perturbations in 
the quadratic v x B term. In SH states, on the other hand, one would require Ohm's 
law to be valid for the helical equilibrium. 

The parallel component of Ohm's law gives, in stationary conditions, 

- V • ($B) + IIb • V0 = 77j • B (57) 

ZTT 

where $ is the electrostatic potential and Vt is the toroidal loop voltage. Multiplying 
by the Jacobian of the coordinate system in use, and integrating over the two angular 
variables, the term involving the electrostatic potential cancels, and one finds 

^(i?*>=r;(j.B). (58) 

This relationship defines the Ohmic constraint that any stationary equilibrium should 
satisfy. 

The two sides of the Ohmic constraint are plotted in Fig. [5] for a typical 1.5 MA 
SHAx state. As for the power balance described in the preceding section, a fiat Zi,ff 
profile has been assumed. A remarkable discrepancy can be seen, with the first term 
being larger than the second one in the inner portion of the plasma, and smaller in the 
outer one. Such discrepancy could be partially resolved assuming a profile of effective 
charge with a peak in the center of the plasma, instead of the flat profile which has 
been assumed here. Even though a screening of impurity influx has been observed 
in Laser Blow-Off experiments [32], a trapping in the central part of the plasma of 
impurities existing prior to the SHAx state onset might occur. Even if this is the 
case, the a — Qq model assumed for the zeroth-order parallel current density appears 
anyway to be not appropriate, since in the outer part of the plasma the profile of {B'^) 
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changes sign, whereas the other one does not. This requires either to assume a profile 
of a changing sign in this region, or a residual dynamo contribution of the secondary 
modes. These considerations point to the need of performing equilibrium calculations 
which take into account the Ohmic constraint. 




Figure 8. Left hand side (open circles) and right hand side (sohd circles) of 
equation II58I I. plotted as a function of the effective radius p. The data refer to 
shot 22182, at t = 49 ms. 



7. Conclusions 

In this paper we have described a perturbative method for the computation of helical 
equilibria for the SHAx states of RFP plasmas in force-free approximation. By using 
the limitation of assuming a zeroth-order axisymmetric equilibrium with circular flux 
surfaces, it has been possible to reduce its calculation, and the subsequent calculation 
of the dominant mode eigenfunction, to the solution of ordinary differential equations. 
The helical flux given by the superposition of the zeroth-order equilibrium and of 
the eigenfunction has been demonstrated to be a good flux function, by comparison 
with spatially resolved measurements of different plasma quantities. Both the force 
free approximation and the perturbative approach do not appear to give errors which 
go beyond the statistical uncertainty of the measurements used for the validation, at 
least for the conditions which have been explored up to now in the RFX-mod device 
(plasma currents up to 1.8 MA). In particular, being the (3 value of these plasmas 
of a few percent, the approximation of neglecting the pressure gradient term in the 
equilibrium equation appears to be valid, despite the presence of a region of strong 
gradient. Large increases in performance are allowed before this approximation breaks 
down. 

As already discussed in ref. H], the possibility of mapping the kinetic 
measurements on the helical flux surfaces resulting from the equilibrium reconstruction 
constitutes a proof of the existence of at least some reminiscence of such surfaces 
in the plasma. While there are no proofs that actual KAM surfaces are formed, 
it has been shown that cantori, that is broken KAM surfaces, are sufficient for 
supporting temperature gradients [33j . These ordered structures have been named 
"ghost surfaces" [34] . The possibility that such ghost surfaces appear in the otherwise 
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chaotic RFX plasma when the SHAx state is attained is the most hkely explanation 
for the experimental observations. 

The availability of an equilibrium reconstruction for the high performance SHAx 
states, as given by the SHEq code which implements the method that we have 
described, opens the path for several lines of research. Among the work in progress 
are the use of the ASTRA code [35] for performing transport calculations taking into 
account the helical geometry and the evaluation of crucial quantities such as the safety 
factor profile. The output of our equilibrium calculation could also serve as the basis 
for stability calculations or for the evaluation of the neoclassical transport coefficients 
in SHAx plasmas. 
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Appendix A 

The covariant metric tensor for a curvilinear coordinate system = (u^^u'^^u^) is 
defined by 

gTj = • e, , (A.l) 

where 

while the contravariant one is defined by 

g^i = Vw* • Vw^ (A.3) 

where 

V.^ = -. (A.4) 

The two tensors are related by g^^ ■ g^j, =5],, so gjf is the inverse matrix of 5"^, . 
The Jacobian of the coordinate system is 



'gZ = A/det[5"] = (Vm^ • Vu^ x Vu^y^ (A.5) 



Given a vector A, its contravariant components are defined as 

A' =A- Vu' or A' = g^^Aj (A.6) 

the second one expresses the contravariant component in terms of the covariant ones. 
In the same way, 

A, =A-e, or A, = A-'' (A.7) 
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Appendix B 

The geometric coordinates = (r, 9, 0) are defined by relationships ([3]) and ([4]). Using 
(jA.ip and (jA.2[) the gfj elements, that link them to the cartesian coordinates, can be 
explicitly calculated: 

/ 1 - 2 A' cos 61 + A'2 rA'sine \ 
= rA'sin6' (B.8) 

\ J 

The gl^ elements of the inverse matrix are usually calculated from the formulas for 
the inversion of a block diagonal matrix (g"^ = = and — g^^^ = 0). The 

Jacobian y/gZ = ^det[gfj] — (Vr • V0 x V(/>)~^ of the geometric coordinates is given 

by 

Vff;:^ ri?(l - A'cos6') (B.9) 

In the text we have defined the coordinate system = (r, i?, 0), which is the fiux 
coordinate system of the zeroth-order axisymmetric equilibrium Bq, as the coordinate 
system that one obtains deforming the poloidal angle 6 of the geometrical coordinates 
in order to achieve straight magnetic field lines: 

^ = + X{r,e) (B.IO) 

From Ampere's law we have also obtained an explicit expression for A(r, 9), (flUl) . The 
contravariant metric tensor elements can be computed using the relations between the 
gradients of the two coordinate systems: 

Vr = Vr (B.ll) 



dr 

V(/> = V(/) (B.13) 



It is therefore found 



while all the other elements are equal to those of the geometric coordinates. The 
g™- elements of the inverse matrix can again be calculated from the formulas for the 
inversion of a block diagonal matrix [g"^ = g'^ = and g'fj' — g'jf = 0), or can be 
computed by writing the relation between the w'^ and the {R, 0, Z) system. In this 
second case, using equations p3)) and (fT4|) the covariant metric tensor elements are 
found to be 

ff^; = 1 + f2A'2 + + ^A"2 - ^A" - |1a' + rA'A") - 2A'cosi? + o«^.16) 
V 2 /to /to / 

9Z = (l + ^ + - ]^^') - - ^^^^ + «(^'^') 

g";;^ = r (^r A" + A' - -^^ sin i9 + oCe^fe) (B.18) 

ffL = i?'- (B.19) 
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'^ith g^^ — g"^^ — g^^ = g^^ = 0. In computing these elements the approximation used 
in ref.[3^ has been adopted of retaining the secular terms (i.e. those not dependent 
on d) up to o(e^) and the harmonics up to o(e). 
The Jacobian is 

1 K{r) 



/9u 



i?2 



(B.20) 



with K{r) already written in (jl6p . In order to compute the harmonics of the first 
order perturbation, one wishes to expand also the metric tensor elements: 

111,., , , 

(B.21) 



:e + c.c. 



with 



/gwo 



'9wi 



= IH J A' + o(e^) 

ffiL.o ^Rq \ Rq Rq 2i?o 

1 1 



/gwi 



(B.22) 



(B.23) 



being the to = O/n = and to = ±l/n = harmonics of the Jacobian respectively. 
We can also relate the Jacobian ^/g^ to J^- Using (jB.12[) one obtains 



1 



1 



1 



(B.24) 



Finally, the metric tensor element combinations appearing in eqs. p7p and psp 



are 



Q 



9m 



1 



l + A'M- + 2r2 



Rq 



4I!. _ A 



'9u, 

^/9^ 

w 
tfrr 

9m 



0,0 



A'2 
"2 



0,0 



= — ( 1 + — + — + —A' - — + o{t^) 

Rq \ 2i?Q 2 2i?o ^ 

o{e%) 



A 



±1,0 



±1,0 



±1,0 



/gwo 
1 

^/9wo \Ro 
1 



A' + o(e3) 



Die 



± 



rA" + A'-^+o{e') 



■km) 

(B.26) 
(B.27) 
(B.28) 
(B.29) 
(B.30) 



Appendix C 

The definition of the helical coordinates (x, used in this paper to describe SHAx 
states is given in the text: the radial coordinate is the helical flux x (HH), the poloidal- 
like angle /3 is given in and the toroidal angle is the usual one. 

The contravariant metric tensor elements are derived from the relation between 
the gradients in this coordinate system and in the flux coordinate system of the zeroth- 
order axisymmetric magnetic field, = (r, i^, 0). By definition the /3 angle is a 
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function of the cylindrical {R, (/>, Z) coordinates, where the toroidal dependence is due 
to the dependence on the toroidal angle of the magnetic axis position, Ra{<t>) and 
Za{(t)), SO that 

Vx = ^Vr + ^Vt? + (C.l) 







dx 


or 




' d4> 






d/3 


or 






Vcj) 







V;3 = ^Vr + ^Vi9 + (C.2) 

(C.3) 

These relationships allow to compute all the components of the contravariant metric 
tensor. Here we write explicitly just the g^^ element (used in the text for the 
calculation of the thermal conductivity): 

dr J \dd J j 9w + J J 9w 

The partial derivatives of x with respect to d and cf) are straightforward to compute, 
whereas the one with respect to r is evaluated through expressions ((37)) and (|38l) . 
The Jacobian is also derived from (|C.3p . following the definition 

V^-(Vx-V/3x V0)"\ (C.4) 

so that 

J=.v.-,v,.v«.(|g-|f)-i=. 

The proof of this simple relation goes as follows: 

dB dB 
or ov 



so 



Vx • V/3 X V(/) = 

^ f ^Vr + ?|V^ + ^ Vp/^z") . Vr X V0 + ^V^ x V0 



9r dcf) J \ dr d-d 

^^Vr • Vi? X V0 + ^^Vi? • Vr X V0 
or ay ov or 

dxdl3 dxdl3\ 

^ Vr • V-d X V(/) 



9r d-d d-d dr ^ 

that is exactly relation (jC.5|) . The partial derivatives of /3 are computed using the 
definition of /3 in terms of R and Z and relations ([T^ and ([H]). 
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